Setup

Load packages

library(here)
library(ggplot2)
library(scales)
library(ggh4x)
library(patchwork)
library(ggchromatic)
library(ggforce)
library(gganimate)

library(DESeq2)
library(SummarizedExperiment)
library(seriation)

Load data

exp <- readRDS(here("rds", "experiment_timeseries.rds"))

Setup Aesthetics

mycolour <- "#000000FF" # Opaque Black

theme_set(theme_gray())
theme_update(text = element_text(colour = mycolour),
             line = element_line(colour = mycolour),
             aspect.ratio = 1,
             axis.line  = element_line(colour = mycolour),
             axis.ticks = element_line(colour = mycolour),
             axis.text  = element_text(colour = mycolour),
             legend.key = element_blank(),
             panel.background = element_blank(),
             panel.grid.major = element_line(colour = "grey95"),
             panel.grid.minor = element_blank(),
             plot.background  = element_blank(),
             strip.background = element_blank(),
             strip.text = element_text(colour = mycolour))

rm(mycolour)

scale_colour_dtag <- function(...) {
  scale_colour_manual(
    name = "dTAG",
    values = c(
      "SOX2"  = "#e41a1c",
      "OCT4"  = "#377eb8",
      "NANOG" = "#4daf4a"
    ),
    ...
  )
}

scale_fill_dtag <- function(...) {
  scale_fill_manual(
    name = "dTAG",
    values = c(
      "SOX2"  = "#e41a1c",
      "OCT4"  = "#377eb8",
      "NANOG" = "#4daf4a"
    ),
    ...
  )
}

fill_shape_guide <- function(...) {
  guide_legend(override.aes = list(shape = 21), ...)
}

Function

make_cat <- function(data, padj_thres = 0.05, larger = "up", smaller= "down", 
                     ns = "n.s.") {
    with(as.data.frame(data), ifelse(is.na(padj) | padj >= padj_thres, ns,
                                   ifelse(log2FoldChange > 0, larger, smaller)))
}

make_contrast <- function(item, exp) {
  as.numeric(resultsNames(exp) %in% item)
}

wrap_ma_plots <- function(result_list) {
  df <- lapply(result_list, as.data.frame)
  df <- data.table::rbindlist(df, idcol = "time")
  df <- transform(
    df,
    cat  = make_cat(df),
    time = factor(time, levels = times)
  )
  
  lims <- range(df$log2FoldChange)
  lims <- c(-1, 1) * max(abs(lims))
  
  g <- ggplot(df, aes(baseMean, log2FoldChange)) +
    geom_point(aes(colour = cat),
               size = 0.1, alpha = 0.1, shape = 16) +
    geom_text(
      aes(
        label = after_stat(number(count)),
        group = cat, colour = cat,
        x = stage(1, after_stat = Inf),
        y = after_stat(rep_len(c(-3, 0, 3), length(x)))
      ),
      stat = "count", inherit.aes = FALSE, hjust = 1,
      show.legend = FALSE
    ) +
    scale_x_continuous(
      name = expression(mu*" Normalised Expression"),
      trans = "log10",
      guide = guide_axis_logticks(),
      limits = c(1, 10000),
      expand = c(0,0),
      labels = number_format(accuracy = 1)
    ) +
    scale_y_continuous(
      limits = c(-1, 1) * max(abs(lims)),
      breaks = breaks_width(2),
      name = expression("Log"[2]*" Fold Change vs 0H")
    ) +
    scale_colour_manual(
      name = "Direction",
      values = c("down" = "#ff5c49", "n.s." = "#949494", "up" = "#009bef"),
      guide = guide_legend(override.aes = list(size = 2, alpha = 1))
    ) +
    facet_wrap(~ time) +
    theme(axis.ticks.length.x = unit(5.5, "pt"),
          aspect.ratio = 2/(1 + sqrt(5)))
  
  list(plot = g, data = df)
}

parallel_sets_plot <- function(data) {
  data$id <- unname(unlist(lapply(split(data$time, data$time), seq_along)))
  ids <- with(data, sort(unique(id[padj < 0.05 & !is.na(padj)])))
  data <- subset(data, id %in% ids)
  
  ggplot(data, aes(time, id = id, split = cat, value = 1)) +
  geom_parallel_sets(
    axis.width = 0.15, alpha = 0.5,
    aes(fill = after_stat(split))
  ) +
  geom_parallel_sets_axes(
    axis.width = 0.15,
    aes(fill = after_stat(label))
  ) +
  geom_parallel_sets_labels(
    aes(label = after_stat(number(value, accuracy = 1)))) +
  scale_fill_manual(
    name = "Direction",
    values = c("down" = "#ff5c49", "n.s." = "#949494", "up" = "#009bef"),
    guide = guide_legend(override.aes = list(size = 2, alpha = 1))
  ) +
  scale_x_discrete(
    name = "Time", expand = c(0, 0.05)
  ) +
  scale_y_continuous(expand = c(0,0)) +
  coord_cartesian(clip = "off") +
  guides(y = "none", x = guide_axis(position = "top")) +
  theme(aspect.ratio = 2/(1 + sqrt(5)),
        panel.grid.major = element_blank())
}

Aim

Do differential accessibility analysis for the timeseries data.

exp$run <- factor(exp$run)
exp <- DESeqDataSet(exp, ~ run + tag + time + time:tag)
## renaming the first element in assays to 'counts'
## converting counts to integer mode

Filter

thres_count <- 10
thres_num <- 6

seqinfo <- seqinfo(rowRanges(exp))
seqinfo <- keepStandardChromosomes(seqinfo, "Mus_musculus")
seqnames <- seqnames(seqinfo)
seqnames <- seqnames[seqnames != "chrM"]

rs <- rowSums(assay(exp) > thres_count)
row_keep <- rs > thres_num
row_keep <- row_keep & seqnames(rowRanges(exp)) %in% seqnames
row_keep <- row_keep & rowMeans(assay(exp)) < 1000

col_keep <- exp$basename != "11_NANOG_0H_r1_TAAGGCGA-TCTTACGC_S1"
# col_keep <- col_keep & exp$tag != "NANOG"

exp <- exp[row_keep, col_keep]

Size factors

exp$tag <- droplevels(exp$tag)
exp <- estimateSizeFactors(exp, type = "poscounts")


df <- data.frame(
  tag = exp$tag,
  time = exp$time,
  rip = colSums(assay(exp)),
  sizefactor = sizeFactors(exp)
)

fit <- MASS::rlm(sizefactor ~ 0 + rip, data = df)

ggplot(df, aes(rip, sizefactor)) +
  geom_abline(slope = coef(fit), colour = "grey50", linetype = 2) +
  geom_point(aes(fill = tag, shape = time),
             size = 3) +
  scale_fill_dtag(guide = fill_shape_guide()) +
  scale_shape_manual(values = c(25, 21:24),
                     name = "Time") +
  scale_x_continuous(
    name = "Reads in Peaks",
    labels = number_format(1, scale = 1e-6, suffix = "M"),
    limits = c(0, NA)
  ) +
  scale_y_continuous(
    name = "Size Factor",
    limits = c(0, NA)
  )

Mean-dispersion trend

Something seems off here. I would have expected the points to be closer to the line, which likely means that there is some unknown variation in the model that isn’t been accounted for.

exp <- estimateDispersions(exp, quiet = TRUE)

df <- as.data.frame(rowData(exp))
df <- df[order(df$baseMean),]

ggplot(df, aes(baseMean, dispMAP)) +
  geom_point(shape = 16, alpha = 0.1, size = 0.1) +
  geom_line(aes(y = dispFit),
            colour = "red", linetype =  2) +
  scale_x_continuous(
    trans = "log10",
    guide = guide_axis_logticks(),
    labels = number_format(1),
    name = expression(mu * " Normalised Accessibility"),
    limits = c(1, 10000),
    expand = c(0,0)
  ) +
  scale_y_continuous(
    trans = "log10",
    guide = guide_axis_logticks(),
    expand = c(0, 0),
    name = expression("Dispersion Estimate "*varphi),
    limits = c(1e-3, 1)
  ) +
  theme(axis.ticks.length = unit(5.5, "pt"))

PCA

norm <- assay(vst(exp))
pca <- prcomp(t(norm))

df <- data.frame(
  pc1 = pca$x[, 1],
  pc2 = pca$x[, 2],
  factor = exp$tag,
  time = exp$time
)

varexp <- with(pca, sdev^2 / sum(sdev^2))

pcplot <- ggplot(df, aes(pc1, pc2)) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_point(aes(fill = factor, shape = time), size = 3) +
  scale_fill_dtag(guide = fill_shape_guide()) +
  scale_shape_manual(values = c(25, 21:24),
                     name = "Time") +
  labs(
    x = paste0("PC1 (", percent(varexp[1], 0.1), ")"),
    y = paste0("PC2 (", percent(varexp[2], 0.1), ")")
  ) +
  coord_equal() +
  theme(aspect.ratio = NULL,
        axis.line = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "bottom",
        legend.box = "vertical",
        legend.box.just = "left",
        legend.box.spacing = unit(0, "cm"),
        legend.box.margin = margin(),
        legend.spacing.y = unit(0, "cm"),
        legend.margin = margin())

vardf <- data.frame(x = seq_along(varexp), y = varexp)

varplot <- ggplot(vardf, aes(x, y)) +
  geom_pointpath() +
  scale_y_continuous(
    limits = c(0, NA),
    expand = c(0,0,0.05,0),
    labels = percent_format(),
    name = "Variance Explained"
  ) +
  scale_x_continuous(limits = c(0, NA),
                     expand = c(0,0),
                     name = "PC #") +
  theme(aspect.ratio = (1 + sqrt(5))/2,
        panel.grid.major.x = element_blank())

pcplot | varplot

Clustering

dist <- as.dist(1 - cor(norm, method = "spearman"))

clust <- hclust(dist, method = "ward.D2")
clust <- reorder(clust, dist = dist, method = "OLO")

df <- reshape2::melt(cor(norm, method = 'spearman'))

ggplot(df, aes(Var1, Var2, fill = value)) +
  geom_raster() +
  guides(y.sec = "axis", x.sec = guide_axis(position = "bottom")) +
  GENOVA:::scale_fill_GENOVA_div(name = expression("Spearman's "*rho)) +
  scale_x_dendrogram(hclust = clust, name = NULL,
                     guide = guide_dendro(label = FALSE, position = "top")) +
  scale_y_dendrogram(hclust = clust, name = NULL,
                     guide = guide_dendro(label = FALSE)) +
  theme(axis.text.x.bottom = element_text(angle = 90, hjust = 1, vjust = 0.5))

Testing differential accessibility

We have setup DESeq2’s design argument as follows:

\[X_{ij} = \alpha_i + \beta_{1,i}\text{run}_j + \beta_{2,i}\text{tag}_j + \beta_{3,i}\text{time}_j + \beta_{4,i}(\text{tag}_j \times \text{time}_j) + \varepsilon\] Wherein:

  • \(X_{ij}\) is the accessibility of peak \(i\) in sample \(j\).
  • \(a_i\) is the intercept term, which will catch the average accessibility of peak \(i\).
  • \(\beta_{1,i}\) is a term for peak \(i\) how much the sequencing run of sample \(j\) will influence peak accessibility \(X_{ij}\). Including this term serves as a sponge for batch effects so that they don’t seep into coefficients we actually care about.
  • \(\text{run}_j\) is a factor variable for the sample run with run number 5831 as reference level and run 5859 and 5980 as binary coefficients.
  • \(\beta_{2,i}\) is a coefficient for peak \(i\) how much the accessibility changes between different dTAG lines at 0H (because 0H is the reference level for \(\text{time}_j\)). The actual coefficient is not really interesting as it captures mostly clonal and parental differences between cell lines.
  • \(\text{tag}_j\) is a factor variable encoding the dTAG line. We change the reference level around for testing differences between lines, whereby the other dTAG lines are encoded as binary coefficients.
  • \(\beta_{3,i}\) is a coefficient for peak \(i\) how much the accessibility changes at each time point in general over all dTAG lines.
  • \(\text{time}_j\) is a factor variable encoding the time point, wherein 0H is the reference level and other timepoints are encoded as binary coefficients.
  • \(\beta_{4,i}\) is an interaction coefficient for how much peak accessibility \(X_{ij}\) changes over time for a specific factor relative to the reference \(\text{tag}_j\) level. This is probably the most interesting coefficient. Note that because \(\text{tag}_j\) and \(\text{time}_j\) are binary coefficents, \(\text{tag}_j \times \text{time}_j\) is 1 if, and only if, sample \(j\) is treated for 2H or longer and the tag is not the reference tag.
  • \(\varexpsilon\) is an error term catching residual variation.

Time

We’re testing the main effects of the time variable using the Wald test with different reference levels (tagged transcription factors).

# We re-level the tag to investigate different dTAGs. This does nothing to the
# stats, but changes what the reference level is. The `time` factor then applies
# for that tag, without having to pick a complex contrast, 
# which simplifies interpretation.
tags <- c(oct4 = "OCT4", sox2 = "SOX2", nanog = "NANOG")

exps <- lapply(tags, function(tag) {
  tmp <- exp
  tmp$tag <- relevel(tmp$tag, tag)
  nbinomWaldTest(tmp)
})

times <- setNames(nm = tail(levels(exp$time), -1))

time <- lapply(exps, function(tmp) {
  lapply(times, function(i) {
    results(tmp, contrast = c("time", i, "0H"))
  })
})

Oct4

oct4 <- wrap_ma_plots(time$oct4)
oct4$plot

parallel_sets_plot(oct4$data)

Sox2

sox2 <- wrap_ma_plots(time$sox2)
sox2$plot

parallel_sets_plot(sox2$data)

Nanog

nanog <- wrap_ma_plots(time$nanog)
nanog$plot

parallel_sets_plot(nanog$data)

Interaction effects

We’ll be doing two tests for interactions. First we do the log-ratio test (LRT), which essentially tests for every gene if it is significantly better modelled by the interaction factor across all time points. We use this as a measure of global interaction.

interact <- list()

Oct4-Sox2

lrt <- exps$oct4
lrt <- lrt[, lrt$tag != "NANOG"]
lrt$tag <- droplevels(lrt$tag)
lrt <- nbinomLRT(lrt, reduced = ~ run + tag + time)
## found results columns, replacing these
interact$oct4sox2_lrt <- results(lrt, contrast = make_contrast("tagSOX2.time2H", lrt))
interact$oct4sox2_2h <- results(exps$oct4, contrast = make_contrast("tagSOX2.time2H", exps$oct4))

df <- data.frame(
  oct4 = time$oct4$`2H`$log2FoldChange,
  sox2 = time$sox2$`2H`$log2FoldChange
)
df <- rbind(df, df)
df$cat <- c(make_cat(interact$oct4sox2_lrt), make_cat(interact$oct4sox2_2h))
df$type <- rep(c("Log Ratio Test", "Wald Test at 2H"), each = nrow(exps$oct4))

lim <- range(c(df[[1]], df[[2]]))

ggplot(df, aes(sox2, oct4, colour = cat)) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_point(size = 0.3, alpha = 0.3, shape = 16) +
  geom_abline(slope = 1, intercept = 0, linetype = 2) +
  geom_label(
    aes(colour = cat,
        label = after_stat(number(count)),
        group = cat,
        x = after_stat(c(-3, 0, 3, -3, 0, 3)),
        y = stage(0, after_stat = c(3, 0, -3, 3, 0, -3))),
    stat = "count", label.size = 0, label.padding = unit(0.1, "lines"),
    inherit.aes = FALSE, show.legend = FALSE
  ) +
  scale_colour_manual(
    values = c(CYRUP::CYRUP(2), "grey50", CYRUP::CYRUP(3)),
    name = "Effect",
    guide = guide_legend(override.aes = list(alpha = 1, size = 2))
  ) +
  scale_y_continuous(
    limits = lim, 
    name = expression(atop("Oct4 2H vs 0H", "Log"[2]*" Fold Change")),
    breaks = breaks_width(2)
  ) +
  scale_x_continuous(
    limits = lim, 
    name = expression(atop("Sox2 2H vs 0H", "Log"[2]*" Fold Change")),
    breaks = breaks_width(2)
  ) +
  facet_wrap(~ type) +
  theme(
    axis.line = element_blank(),
    axis.ticks = element_blank()
  )

Oct4-Nanog

lrt <- exps$oct4
lrt <- lrt[, lrt$tag != "SOX2"]
lrt$tag <- droplevels(lrt$tag)
lrt <- nbinomLRT(lrt, reduced = ~ run + tag + time)
## found results columns, replacing these
interact$oct4nanog_lrt <- results(lrt, contrast = make_contrast("tagNANOG.time2H", lrt))
interact$oct4nanog_2h <- results(exps$oct4, contrast = make_contrast("tagNANOG.time2H", exps$oct4))

df <- data.frame(
  oct4 = time$oct4$`2H`$log2FoldChange,
  nanog = time$nanog$`2H`$log2FoldChange
)
df <- rbind(df, df)
df$cat <- c(make_cat(interact$oct4nanog_lrt), make_cat(interact$oct4nanog_2h))
df$type <- rep(c("Log Ratio Test", "Wald Test at 2H"), each = nrow(exps$oct4))

lim <- range(c(df[[1]], df[[2]]))

ggplot(df, aes(nanog, oct4, colour = cat)) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_point(size = 0.3, alpha = 0.3, shape = 16) +
  geom_abline(slope = 1, intercept = 0, linetype = 2) +
  geom_label(
    aes(colour = cat,
        label = after_stat(number(count)),
        group = cat,
        x = after_stat(c(-3, 0, 3, -3, 0, 3)),
        y = stage(0, after_stat = c(3, 0, -3, 3, 0, -3))),
    stat = "count", label.size = 0, label.padding = unit(0.1, "lines"),
    inherit.aes = FALSE, show.legend = FALSE
  ) +
  scale_colour_manual(
    values = c(CYRUP::CYRUP(2), "grey50", CYRUP::CYRUP(3)),
    name = "Effect",
    guide = guide_legend(override.aes = list(alpha = 1, size = 2))
  ) +
  scale_y_continuous(
    limits = lim, 
    name = expression(atop("Oct4 2H vs 0H", "Log"[2]*" Fold Change")),
    breaks = breaks_width(2)
  ) +
  scale_x_continuous(
    limits = lim, 
    name = expression(atop("Nanog 2H vs 0H", "Log"[2]*" Fold Change")),
    breaks = breaks_width(2)
  ) +
  facet_wrap(~ type) +
  theme(
    axis.line = element_blank(),
    axis.ticks = element_blank()
  )

Sox2-Nanog

lrt <- exps$sox2
lrt <- lrt[, lrt$tag != "OCT4"]
lrt$tag <- droplevels(lrt$tag)
lrt <- nbinomLRT(lrt, reduced = ~ run + tag + time)
## found results columns, replacing these
interact$sox2nanog_lrt <- results(lrt, contrast = make_contrast("tagNANOG.time2H", lrt))
interact$sox2nanog_2h <- results(exps$sox2, contrast = make_contrast("tagNANOG.time2H", exps$sox2))

df <- data.frame(
  sox2 = time$sox2$`2H`$log2FoldChange,
  nanog = time$nanog$`2H`$log2FoldChange
)
df <- rbind(df, df)
df$cat <- c(make_cat(interact$sox2nanog_lrt), make_cat(interact$sox2nanog_2h))
df$type <- rep(c("Log Ratio Test", "Wald Test at 2H"), each = nrow(exps$sox2))

lim <- range(c(df[[1]], df[[2]]))

ggplot(df, aes(nanog, sox2, colour = cat)) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_point(size = 0.3, alpha = 0.3, shape = 16) +
  geom_abline(slope = 1, intercept = 0, linetype = 2) +
  geom_label(
    aes(colour = cat,
        label = after_stat(number(count)),
        group = cat,
        x = after_stat(c(-3, 0, 3, -3, 0, 3)),
        y = stage(0, after_stat = c(3, 0, -3, 3, 0, -3))),
    stat = "count", label.size = 0, label.padding = unit(0.1, "lines"),
    inherit.aes = FALSE, show.legend = FALSE
  ) +
  scale_colour_manual(
    values = c(CYRUP::CYRUP(2), "grey50", CYRUP::CYRUP(3)),
    name = "Effect",
    guide = guide_legend(override.aes = list(alpha = 1, size = 2))
  ) +
  scale_y_continuous(
    limits = lim, 
    name = expression(atop("Sox2 2H vs 0H", "Log"[2]*" Fold Change")),
    breaks = breaks_width(2)
  ) +
  scale_x_continuous(
    limits = lim, 
    name = expression(atop("Nanog 2H vs 0H", "Log"[2]*" Fold Change")),
    breaks = breaks_width(2)
  ) +
  facet_wrap(~ type) +
  theme(
    axis.line = element_blank(),
    axis.ticks = element_blank()
  )

Notes

pthres <- function(x, thres = 0.05) {!is.na(x) & x < thres}

df <- data.frame(
  oct_lfc = time$oct4$`2H`$log2FoldChange,
  oct_pad = time$oct4$`2H`$padj,
  sox_lfc = time$sox2$`2H`$log2FoldChange,
  sox_pad = time$sox2$`2H`$log2FoldChange,
  int_lfc = interact$oct4sox2_2h$log2FoldChange,
  int_pad = interact$oct4sox2_2h$padj
)

df$cat <- elzo_cat <- with(df, dplyr::case_when(
  sox_lfc < -1 & oct_lfc < -1 & pthres(oct_pad) & pthres(sox_pad) & !pthres(int_pad) ~ "both_down",
  int_lfc > 0 & pthres(int_pad) & sox_lfc > -0.5 ~ "oct_down",
  int_lfc < 0 & pthres(int_pad) & oct_lfc > -0.5 ~ "sox_down",
  TRUE ~ "n.s."
))

ggplot(df, aes(sox_lfc, oct_lfc, colour = cat)) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_point(shape = 16, alpha = 0.5, size = 0.5) +
  annotate("line", x = c(-Inf, -0.75, -0.75), y = c(-0.75, -0.75, -Inf),
           linetype = 2) +
  annotate("line", x = c(-0.75, Inf), y = c(-0.75, Inf), linetype = 2) +
  geom_label(aes(colour = cat,
                 label = after_stat(count),
                 group = cat,
                 x = stage(0, after_stat = c(-1, 1, 1, -1)*4),
                 y = after_stat(c(-1, 1, -1, 1)*4)),
             stat = "count", label.size = 0, label.padding = unit(0.1, "lines"),
             inherit.aes = FALSE, show.legend = FALSE) +
  scale_colour_manual(
    values = c(CYRUP::CYRUP(1:3), "grey50")[c(1,4,2,3)],
    name = "Category",
    guide = guide_legend(override.aes = list(size = 2, alpha = 1)),
    labels = c("Both", "Neither", "Oct4", "Sox2")
  ) +
  scale_x_continuous(
    breaks = breaks_width(2),
    name = expression(atop("Sox2 2H vs 0H","Log"[2]*" Fold Change"))
  ) +
  scale_y_continuous(
    breaks = breaks_width(2),
    name = expression(atop("Oct4 2H vs 0H", "Log"[2]*" Fold Change"))
  ) +
  theme(axis.line = element_blank(),
        axis.ticks = element_blank(),
        legend.key.size = unit(1, "lines"),
        legend.background = element_blank())

i <- seq_len(nrow(exps$oct4))
get_lfc <- function(x) {
  as.vector(vapply(x, function(y) y$log2FoldChange[i], as.numeric(i)))
}



df <- data.frame(
  oct_lfc = get_lfc(time$oct4),
  sox_lfc = get_lfc(time$sox2),
  time = rep(factor(c("2H", "4H", "6H", "24H"),
                    levels = c("2H", "4H", "6H", "24H")),
             each = length(i)),
  cat = rep(elzo_cat[i], 4),
  grp = rep(seq_along(i), 4)
)

# ii <- sample(i, 100)
# df <- df[df$grp %in% ii,]

lim <- range(c(df$oct_lfc, df$sox_lfc))

lines <- data.frame(
  x = c(-Inf, -0.75, -0.75, -0.75, Inf),
  y = c(-0.75, -0.75, -Inf, -0.75, Inf),
  grp = c(1, 1, 1, 2, 2)
)

plt <- ggplot(df, aes(sox_lfc, oct_lfc)) +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  geom_point(
    aes(colour = cat, group = grp),
    shape = 16, size = 0.5, alpha = 0.5
  ) +
  geom_path(
    data = lines,
    aes(x, y, group = grp),
    linetype = 2
  ) +
  scale_colour_manual(
    values = c(CYRUP::CYRUP(1:3), "grey50")[c(1,4,2,3)],
    name = "Category",
    guide = guide_legend(override.aes = list(size = 2, alpha = 1)),
    labels = c("Both", "Neither", "Oct4", "Sox2")
  ) +
  scale_x_continuous(
    limits = lim, breaks = breaks_width(2),
    name = expression(atop("Sox2 2H vs 0H", "Log"[2]*" Fold Change"))
  ) +
  scale_y_continuous(
    limits = lim, breaks = breaks_width(2),
    name = expression(atop("Oct4 2H vs 0H", "Log"[2]*" Fold Change"))
  ) +
  theme(axis.line = element_blank(),
        axis.ticks = element_blank(),
        legend.key.size = unit(1, "lines"),
        legend.background = element_blank())

anim <- plt + transition_states(time, state_length = 0.01, wrap = FALSE) +
  ease_aes("cubic-in-out") +
  ggtitle("At {closest_state}") +
  shadow_wake(0.2, wrap = FALSE)

anim <- animate(
  anim, renderer = gifski_renderer(),
  device = "ragg_png", res = 150, width = 1000, height = 1000
)

anim_save(here("export", "da_anim.gif"), animation = anim)
knitr::include_graphics(here("export", "da_anim.gif"))

results <- c(
  list(rowRanges = rowRanges(exps$oct4)),
  time,
  interact
)

saveRDS(results, here("rds", "01_DA_Results.rds"))

References

Session Info

Platform
settings values
version R version 3.6.3 (2020-02-29)
os Ubuntu 16.04.7 LTS
system x86_64, linux-gnu
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz Europe/Amsterdam
date 2021-04-06
user t.vd.brand
Packages
package attached version date source
Attached base packages
base ★ 3.6.3 2020-02-29 local
compiler 3.6.3 2020-02-29 local
datasets ★ 3.6.3 2020-02-29 local
graphics ★ 3.6.3 2020-02-29 local
grDevices ★ 3.6.3 2020-02-29 local
grid 3.6.3 2020-02-29 local
methods ★ 3.6.3 2020-02-29 local
parallel ★ 3.6.3 2020-02-29 local
splines 3.6.3 2020-02-29 local
stats ★ 3.6.3 2020-02-29 local
stats4 ★ 3.6.3 2020-02-29 local
tools 3.6.3 2020-02-29 local
utils ★ 3.6.3 2020-02-29 local
Attached other packages
Biobase ★ 2.46.0 2019-10-29 Bioconductor
BiocGenerics ★ 0.32.0 2019-10-29 Bioconductor
BiocParallel ★ 1.20.1 2019-12-21 Bioconductor
DelayedArray ★ 0.12.3 2020-04-09 Bioconductor
DESeq2 ★ 1.26.0 2019-10-29 Bioconductor
GenomeInfoDb ★ 1.22.1 2020-03-27 Bioconductor
GenomicRanges ★ 1.38.0 2019-10-29 Bioconductor
gganimate ★ 1.0.7 2020-10-15 CRAN (R 3.6.3)
ggchromatic ★ 0.0.1.0 2021-02-24 local
ggforce ★ 0.3.2 2020-06-23 CRAN (R 3.6.3)
ggh4x ★ 0.1.2.1.9 2021-03-31 local
ggplot2 ★ 3.3.2.9000 2020-07-11 local
here ★ 1.0.0 2020-11-15 CRAN (R 3.6.3)
IRanges ★ 2.20.2 2020-01-13 Bioconductor
matrixStats ★ 0.57.0 2020-09-25 CRAN (R 3.6.3)
patchwork ★ 1.1.0 2020-11-09 CRAN (R 3.6.3)
S4Vectors ★ 0.24.4 2020-04-09 Bioconductor
scales ★ 1.1.1 2020-05-11 CRAN (R 3.6.3)
seriation ★ 1.2-9 2020-10-01 CRAN (R 3.6.3)
SummarizedExperiment ★ 1.16.1 2019-12-19 Bioconductor
Loaded via namespace only
annotate 1.64.0 2019-10-29 Bioconductor
AnnotationDbi 1.48.0 2019-10-29 Bioconductor
backports 1.2.0 2020-11-02 CRAN (R 3.6.3)
base64enc 0.1-3 2015-07-28 CRAN (R 3.6.1)
Biostrings 2.54.0 2019-10-29 Bioconductor
bit 4.0.4 2020-08-04 CRAN (R 3.6.3)
bit64 4.0.5 2020-08-30 CRAN (R 3.6.3)
bitops 1.0-6 2013-08-17 CRAN (R 3.6.1)
blob 1.2.1 2020-01-20 CRAN (R 3.6.2)
checkmate 2.0.0 2020-02-06 CRAN (R 3.6.2)
cluster 2.1.1 2021-02-14 CRAN (R 3.6.3)
codetools 0.2-18 2020-11-04 CRAN (R 3.6.3)
colorspace 2.0-0 2020-11-11 CRAN (R 3.6.3)
crayon 1.3.4 2017-09-16 CRAN (R 3.6.1)
CYRUP 0.9 2019-07-02 Github ()
data.table 1.13.4 2020-12-08 CRAN (R 3.6.3)
DBI 1.1.0 2019-12-15 CRAN (R 3.6.1)
digest 0.6.27 2020-10-24 CRAN (R 3.6.3)
dplyr 1.0.2 2020-08-18 CRAN (R 3.6.3)
ellipsis 0.3.1 2020-05-15 CRAN (R 3.6.3)
evaluate 0.14 2019-05-28 CRAN (R 3.6.1)
farver 2.0.3.9000 2021-01-30 Github ()
foreach 1.5.1 2020-10-15 CRAN (R 3.6.3)
foreign 0.8-76 2020-03-03 CRAN (R 3.6.3)
Formula 1.2-4 2020-10-16 CRAN (R 3.6.3)
genefilter 1.68.0 2019-10-29 Bioconductor
geneplotter 1.64.0 2019-10-29 Bioconductor
generics 0.1.0 2020-10-31 CRAN (R 3.6.3)
GenomeInfoDbData 1.2.2 2019-11-19 Bioconductor
GenomicAlignments 1.22.1 2019-11-12 Bioconductor
GENOVA 1.0.0 2021-02-10 local
ggdendro 0.1.22 2020-09-13 CRAN (R 3.6.3)
gifski 1.4.3 2021-03-22 CRAN (R 3.6.3)
glue 1.4.2 2020-08-27 CRAN (R 3.6.3)
gridExtra 2.3 2017-09-09 CRAN (R 3.6.1)
gtable 0.3.0 2019-03-25 CRAN (R 3.6.1)
Hmisc 4.4-1 2020-08-10 CRAN (R 3.6.3)
hms 0.5.3 2020-01-08 CRAN (R 3.6.2)
htmlTable 2.1.0 2020-09-16 CRAN (R 3.6.3)
htmltools 0.5.0 2020-06-16 CRAN (R 3.6.3)
htmlwidgets 1.5.2 2020-10-03 CRAN (R 3.6.3)
iterators 1.0.13 2020-10-15 CRAN (R 3.6.3)
jpeg 0.1-8.1 2019-10-24 CRAN (R 3.6.2)
knitr 1.30 2020-09-22 CRAN (R 3.6.3)
labeling 0.3 2014-08-23 CRAN (R 3.4.0)
lattice 0.20-41 2020-04-02 CRAN (R 3.6.3)
latticeExtra 0.6-29 2019-12-19 CRAN (R 3.6.2)
libTvdB 0.1.0 2020-06-22 local
lifecycle 0.2.0 2020-03-06 CRAN (R 3.6.2)
locfit 1.5-9.4 2020-03-25 CRAN (R 3.6.3)
magrittr 1.5 2014-11-22 CRAN (R 3.4.0)
MASS 7.3-53.1 2021-02-12 CRAN (R 3.6.3)
Matrix 1.3-2 2021-01-06 CRAN (R 3.6.3)
memoise 1.1.0 2017-04-21 CRAN (R 3.6.1)
munsell 0.5.0 2018-06-12 CRAN (R 3.6.1)
nnet 7.3-15 2021-01-24 CRAN (R 3.6.3)
pillar 1.4.7 2020-11-20 CRAN (R 3.6.3)
pkgconfig 2.0.3 2019-09-22 CRAN (R 3.6.1)
plyr 1.8.6 2020-03-03 CRAN (R 3.6.2)
png 0.1-7 2013-12-03 CRAN (R 3.6.2)
polyclip 1.10-0 2019-03-14 CRAN (R 3.6.1)
prettyunits 1.1.1 2020-01-24 CRAN (R 3.6.2)
progress 1.2.2 2019-05-16 CRAN (R 3.6.1)
purrr 0.3.4 2020-04-17 CRAN (R 3.6.3)
R6 2.5.0 2020-10-28 CRAN (R 3.6.3)
RColorBrewer 1.1-2 2014-12-07 CRAN (R 3.4.0)
Rcpp 1.0.5 2020-07-06 CRAN (R 3.6.3)
RCurl 1.98-1.2 2020-04-18 CRAN (R 3.6.3)
registry 0.5-1 2019-03-05 CRAN (R 3.6.1)
reshape2 1.4.4 2020-04-09 CRAN (R 3.6.3)
rlang 0.4.8 2020-10-08 CRAN (R 3.6.3)
rmarkdown 2.6 2020-12-14 CRAN (R 3.6.3)
rpart 4.1-15 2019-04-12 CRAN (R 3.6.1)
rprojroot 2.0.2 2020-11-15 CRAN (R 3.6.3)
Rsamtools 2.2.3 2020-02-23 Bioconductor
RSQLite 2.2.1 2020-09-30 CRAN (R 3.6.3)
rstudioapi 0.13 2020-11-12 CRAN (R 3.6.3)
rtracklayer 1.46.0 2019-10-29 Bioconductor
stringi 1.5.3 2020-09-09 CRAN (R 3.6.3)
stringr 1.4.0 2019-02-10 CRAN (R 3.6.1)
survival 3.2-7 2020-09-28 CRAN (R 3.6.3)
tibble 3.0.4 2020-10-12 CRAN (R 3.6.3)
tidyselect 1.1.0 2020-05-11 CRAN (R 3.6.3)
TSP 1.1-10 2020-04-17 CRAN (R 3.6.3)
tweenr 1.0.1 2018-12-14 CRAN (R 3.6.1)
vctrs 0.3.5 2020-11-17 CRAN (R 3.6.3)
withr 2.3.0 2020-09-22 CRAN (R 3.6.3)
xfun 0.19 2020-10-30 CRAN (R 3.6.3)
XML 3.99-0.3 2020-01-20 CRAN (R 3.6.2)
xtable 1.8-4 2019-04-21 CRAN (R 3.6.1)
XVector 0.26.0 2019-10-29 Bioconductor
yaml 2.2.1 2020-02-01 CRAN (R 3.6.2)
zlibbioc 1.32.0 2019-10-29 Bioconductor